library(tidyverse)
library(ggplot2)
library(survey)
library(srvyr)
library(pander)
library(knitr)
library(readr)
library(remotes)
library(pewmethods)
library(foreign)
library(haven)
library(estimatr)
library(car)
library(plotrix)
library(xtable)
library(ltm)
library(texreg)
library(reshape2)
library(mice)
library(matrixStats)
library(survey)

options(digits=2)

#You will need to set your own WD
#We suggest dumping the replication files into a new dedicated folder
#And creating an R project set in that folder

#Read data
df_labels<- readRDS("replication_laps.rds")


#############################################################
# Functions 
#############################################################

the.plot.single  <- function(x,ymin=-40,ymax=70,my.axis=T,my.main=NULL,abh0=T,abhcustom=NA,ybsl=F){
  require(plotrix)
  par(mar=c(2,2,1.8,1))
  parties <- length(x$coefficients)
  if(parties==7){xs <- c(4,5,7,6,3,1,2)
  the.pchs <- c(21,25,25,25,24,24,24)
  the.bgs <- c("white",gray(0.7),gray(0.3),1,gray(0.3),gray(0.7),1)}
  if(parties==5){xs <- c(3,4,5,2,1)
  the.pchs <- c(21,25,25,24,24)
  the.bgs <- c("white",gray(0.6),1,gray(0.5),1)}
  if(parties==3){xs <- c(2,3,1) 
  the.pchs <- c(21,25,24)
  the.bgs <- c("white",1,1)}
  plot(c(.5,(parties+0.5)),c(ymin,ymax),type="n"
       ,xlab="",ylab="",xaxt="n",yaxt="n",bty="n")
  mtext(side=3,line=-.5,text=my.main)
  axis(side=2,at=c(ymin,0,ymax),labels=c(ymin,0,ymax))
  if(ybsl==T){axis(side=2,at= x$coefficients[1]
                   ,labels=round(x$coefficients[1],2), tick=T) }
  #abline(h=  x$coefficients[1] ,col=gray(0.7),lty=3)#ploting over the whole area in knit (WTF)
  segments(x0=0.5,x1=(parties+0.5),y0=x$coefficients[1],y1=x$coefficients[1],col=gray(0.7),lty=3)
  if(abh0==T){#abline(h=0,col=gray(0.4),lty=1)
    segments(x0=0.5,x1=(parties+0.5),y0=0,y1=0,col=gray(0.4),)
  }
  #abline(h=abhcustom)#ploting over the whole area in knit (WTF)
  segments(x0=0.5,x1=(parties+0.5),y0=abhcustom,y1=abhcustom)
  plotCI(xs,x$coefficients,
         ui=x$conf.high,li=x$conf.low
         ,col=c(1,gray(0.25),gray(0.25),gray(0.25),1),lwd=1.5
         ,pch=NA,ylab="",xlab="",yaxt="n",xaxt="n",bty="n",add=T)
  points(xs,x$coefficients,pch=the.pchs 
         ,bg=the.bgs
         ,col=1)
}

extract.estse <- function(varname,data=SurveyDesign){
  tmp <- svymean(formula(paste("~",varname)), data, na.rm = T)
  out.est <- as.numeric(tmp[1])
  out.se <- sqrt(attr(tmp,"var"))
  out.ub <- out.est  + 1.96* out.se
  out.lb <- out.est  - 1.96* out.se
  out <- c(est=out.est,se=out.se,ub=out.ub,lb=out.lb)
  return(out)
}


issue.sum <- function(the.issue,invert.polarity=NULL,issue.label=NULL,the.data=df_labels){
  if(is.null(issue.label)){issue.label<-the.issue}
  #Build variable names for the selected issues
  v.agree <- paste(the.issue,"1_agree",sep="")
  v.importance <- paste(the.issue,"2_importance",sep="")
  v.gapantipt <- paste(the.issue,"3_gapantipt",sep="")
  v.gappt <- paste(the.issue,"3_gappt",sep="")
  #Get actual variable: Recoding categorical to binary and keeping numeric
  if(invert.polarity==F){#agreement is conservative
    the.data$y.agree <- the.data[,v.agree]=="Concorda muito"|
      the.data[,v.agree]=="Concorda um pouco"
    the.data$y.gapantipt <-  the.data[[v.gapantipt]]
    the.data$y.gappt <-  the.data[[v.gappt]]
  }else{#disagreement is conservative, use disagreement and get complementary share
    the.data$y.agree <- the.data[,v.agree]=="Discorda um pouco"|
      the.data[,v.agree]=="Discorda muito"    
    the.data$y.gapantipt <-  100-the.data[[v.gapantipt]]
    the.data$y.gappt <-  100-the.data[[v.gappt]]
  }
  #added later... eliminate "nao sabe" from computations of individual agreement
  the.data$y.agree[which(the.data[,v.agree]=="Não sei")] <- NA
  #end addition 
  the.data$y.importance <- the.data[,v.importance]=="Muito importante"
  the.data$y.issue.gap <- the.data$y.gappt-the.data$y.gapantipt  #PT-ANTIPT
  #Check variation across groups
  reg.agree <- lm_robust(y.agree ~type-1,data=the.data,weights=the.data$`weights#all`)
  reg.agree.1 <- lm_robust(y.agree ~type5-1,data=the.data
                           ,weights=the.data$`weights#all`)
  reg.importance <- lm_robust(y.importance ~type5-1,data=the.data
                              ,weights=the.data$`weights#all`)
  reg.importance.1 <- lm_robust(y.importance ~type5-1
                                ,data=the.data,weights=the.data$`weights#all`)
  reg.antipt <- lm_robust(y.gapantipt ~type5-1
                          ,data=the.data,weights=the.data$`weights#all`)
  reg.pt <- lm_robust(y.gappt ~type5-1,data=the.data,weights=the.data$`weights#all`)
  # Gap in preferences: petistas - antipetistas
  require(car)  
  tmp <- linearHypothesis(reg.agree,"typepetista=typeantipetista") #PT-ANTIPT
  real <- as.numeric(attr(tmp,"value"))
  real.se <- sqrt(as.numeric(attr(tmp,"vcov")))
  real.cilow <- real-1.96*real.se
  real.cihigh <- real+1.96*real.se
  real.gap <- c(est=real,se=real.se,cilow=real.cilow,cihigh=real.cihigh)
  the.data$y.relissue.gap <-  log(1+abs(the.data$y.issue.gap) / abs(real.gap["est"]*100))  #THINK HERE
  rm(tmp)
  # Average perception gap
  perceived <- mean(the.data$y.issue.gap,na.rm=T)
  perceived.se <- sd(the.data$y.issue.gap,na.rm=T)
  perceived.cilow <-  perceived-1.96* perceived.se
  perceived.cihigh <-  perceived+1.96* perceived.se
  perceived.gap <- c(est=perceived,se=perceived.se,
                     cilow=perceived.cilow,cihigh=perceived.cihigh)
  #Perception gaps by group
  gap.reg <- lm_robust(y.issue.gap  ~type5-1,data=the.data,
                       weights=the.data$`weights#all`)
  rel.gap.reg <- lm_robust(y.relissue.gap~type5-1,data=the.data
                           ,weights=the.data$`weights#all`)
  #Output
  out <- list(issue.label=issue.label,reg.agree=reg.agree,
              reg.importance=reg.importance,reg.antipt=reg.antipt,
              reg.pt=reg.pt,
              real.gap=real.gap, perceived.gap=perceived.gap,
              perception.gap=gap.reg,
              rel.perception.gap=rel.gap.reg,
              recoded.pt= the.data$y.gappt,
              recoded.antipt=the.data$y.gapantipt)
  return(out)
}

issue.pol <- function(the.issue,invert.polarity=NULL,issue.label=NULL,the.data=df_labels){
  if(is.null(issue.label)){issue.label<-the.issue}
  v.agree <- paste(the.issue,"1_agree",sep="")
  if(invert.polarity==T){
    the.variable <- dplyr::recode(the.data[,v.agree][[1]],`Concorda muito`=-2,`Concorda um pouco`=-1,`Não sei`=0
                                  ,`Discorda um pouco`=1,`Discorda muito`=2
                                  , .default = as.numeric(NA_character_)) 
  }
  if(invert.polarity==F){
    the.variable  <- dplyr::recode(the.data[,v.agree][[1]],`Concorda muito`=2,`Concorda um pouco`=1,`Não sei`=0
                                   ,`Discorda um pouco`=-1,`Discorda muito`=-2
                                   , .default = as.numeric(NA_character_))   
  }
  ## Group polarization 
  the.means <- as.matrix(by(the.variable,df_labels$type,mean,na.rm=T))[c("petista","antipetista"),]
  the.vars <- as.matrix(by(the.variable,df_labels$type,var,na.rm=T))[c("petista","antipetista"),]
  issue.polarization <- (the.means["antipetista"]-the.means["petista"]) / sqrt(the.vars[1]+ the.vars[2])
  
  ## Strong partisan polarization
  the.means.strong <- as.matrix(by(the.variable,df_labels$type5,mean,na.rm=T))[c("petista.strong","antipetista.strong"),]
  the.vars.strong <- as.matrix(by(the.variable,df_labels$type5,var,na.rm=T))[c("petista.strong","antipetista.strong"),]
  issue.polarization.strong <- (the.means.strong["antipetista.strong"]-the.means.strong["petista.strong"]) / 
    sqrt(the.vars.strong["petista.strong"]+ the.vars.strong["antipetista.strong"])
  
  ## Weak partisan polarization
  the.means.weak <- as.matrix(by(the.variable,df_labels$type5,mean,na.rm=T))[c("petista","antipetista"),]
  the.vars.weak <- as.matrix(by(the.variable,df_labels$type5,var,na.rm=T))[c("petista","antipetista"),]
  issue.polarization.weak <- (the.means.weak["antipetista"]-the.means.weak["petista"]) / 
    sqrt(the.vars.weak["petista"]+ the.vars.weak["antipetista"])
  pol <- c(pol=as.numeric(issue.polarization),
           pol.strong=as.numeric(issue.polarization.strong),
           pol.weak=as.numeric(issue.polarization.weak))
}  ## End functions 




#########################################
# Code starts here
# Create partisan groups for comparison
########################################

df_labels$ptID <- car::recode(df_labels$polpartyposit,"'PT'=T;else=F")
df_labels$ptIDstrong <- df_labels$ptID==T & df_labels$petista == "Sim"
df_labels$ptIDextreme <- df_labels$ptID==T & df_labels$petista == "Sim" &
  (df_labels$frequspetista=="Sempre"|df_labels$frequspetista=="Frequentemente") &
  (df_labels$freqpetista=="Sempre"|df_labels$freqpetista=="Frequentemente")

df_labels$antiptID <- car::recode(df_labels$polpartynegat_PT,"'Sim'=T;else=F")
df_labels$antiptIDstrong <- df_labels$antiptID==T & df_labels$antipetista == "Sim"
df_labels$antiptIDextreme <- df_labels$antiptID==T & df_labels$antipetista == "Sim" &
  (df_labels$frequsantipetista=="Sempre"|df_labels$frequsantipetista=="Frequentemente") &
  (df_labels$freqantipetista=="Sempre"|df_labels$freqantipetista=="Frequentemente")

# This is our residual category (includes nonpartisans and partisans of others)
df_labels$noID <- df_labels$ptID == F & df_labels$antiptID == F

# This is an alternative: strict nonpartisans
df_labels$noIDstrict <- df_labels$antiptID==F & is.na(df_labels$polpartyposit)

# Create variables with classification of respondent partisan types. 
# These are used later in several comparisons
df_labels$type <- ifelse(as.numeric(df_labels$ptID)==2,"petista",
                         ifelse(as.numeric(df_labels$antiptID)==2,
                                "antipetista","_nonpartisan"))
df_labels$type5 <- ifelse(as.numeric(df_labels$ptIDstrong),"petista.strong",                                            ifelse(as.numeric(df_labels$antiptIDstrong)                                          ,"antipetista.strong",df_labels$type))
df_labels$type7 <- ifelse(as.numeric(df_labels$ptIDextreme),"petista.extreme",                                            ifelse(as.numeric(df_labels$antiptIDextreme)                                          ,"antipetista.extreme",df_labels$type5))

# Create variables that indicate belonging to certain groups of interest.
# These are used later in several comparisons. 
df_labels$strNE <- df_labels$region=="NE"
df_labels$strCO <- df_labels$region=="CO"
df_labels$str02SM <- df_labels$socialclass =="D-E"
df_labels$str0510SM <- df_labels$socialclass == "B1"| df_labels$socialclass == "B2"
df_labels$stratheist <- df_labels$religionpond =="Ateu, não tem religião"|
  df_labels$religionpond =="É religioso mas não segue nenhuma/ Agnóstico" # the religiao_cotas var includes "prefiro nao responder"in "sem religiao"
df_labels$strcatholic <- df_labels$religiao_cotas =="Católica"
df_labels$strwomen <- df_labels$gender =="Mulher"
df_labels$strblack <- df_labels$colorpond =="Preto"|
  df_labels$colorpond =="Pardo"
df_labels$strwhite <- df_labels$colorpond =="Branco"
df_labels$strSSE <- df_labels$region=="S"|df_labels$region=="SE"
df_labels$str20SM <- df_labels$socialclass =="A"
df_labels$strprotestant <- df_labels$religiao_cotas =="Evangélica"
df_labels$strmen <- df_labels$gender =="Homem"
df_labels$strover60 <-  df_labels$age>=65  
df_labels$str3035 <-  df_labels$age>=30&df_labels$age<=35 
df_labels$catAB <-   df_labels$socialclass == "B1"| df_labels$socialclass == "B2" | df_labels$socialclass == "A"
df_labels$catyoung <- df_labels$age <=35 

# Stratified Independent Sampling design
# Used in some weighted tables and analyses
SurveyDesign <- svydesign(ids = ~1, 
                          weights = ~`weights#all`,
                          strata = NULL,
                          fpc = NULL,
                          fpctype = "correction",
                          data = df_labels)







################################################
# Figure B1 Appendix 
################################################
tp <- get_totals(var = "type5", df=df_labels, wt = "weights#all")[,2]
tp <- rbind(c(tp[5],0,tp[3]),
            tp[c(4,1,2)])
barplot(tp,names.arg=c("Petista","Não Partidário","Antipetista")
        , ylim=c(0,50),ylab=c("Percentage of the Respondents")
        , legend.text = c("Strong sympathizer","Sympathizer")
        , args.legend = list(bty="n",inset=-0.1,xpd=NA))

print(get_totals(var = "type5", df=df_labels, wt = "weights#all"))

#############################################################
# Figure C1  Appendix 
#############################################################

regNE <- lm_robust(strNE~type-1,data=df_labels,weights=df_labels$`weights#all`)
regprotestant <- lm_robust(strprotestant~type-1,data=df_labels,weights=df_labels$`weights#all`)
regwomen  <- lm_robust(strwomen ~type-1,data=df_labels,weights=df_labels$`weights#all`)
regblack <- lm_robust(strblack~type-1,data=df_labels,weights=df_labels$`weights#all`)
regAB <- lm_robust(catAB~type-1,data=df_labels,weights=df_labels$`weights#all`)
regyoung     <- lm_robust(catyoung  ~type-1,data=df_labels,weights=df_labels$`weights#all`)

nr <- 2
nc <- 3
par(oma=c(3,0.1,3,0),mar=c(0,0,0,0),mfrow=c(nr,nc))
the.plot.single(regNE,ymin=0,ymax=.8,my.main="Northeast",abh0=F,ybsl=T)
the.plot.single(regprotestant,ymin=0,ymax=.8,my.main="Protestant",abh0=F,ybsl=T)
the.plot.single(regwomen,ymin=0,ymax=.8,my.main="Women",abh0=F,ybsl=T)
the.plot.single(regblack,ymin=0,ymax=.8,my.main="Non-white",abh0=F,ybsl=T)
the.plot.single(regAB,ymin=0,ymax=.8,my.main="Class AB",abh0=F,ybsl=T)
the.plot.single(regyoung,ymin=0,ymax=.8,my.main="Young",abh0=F,ybsl=T)
title(main = "Types by selected sociodemographic categories",
      outer = TRUE, line = 0.7, cex=2)
par(mfrow=c(1, 1), oma=rep(0, 4), mar=rep(0, 4), new=TRUE)
plot(0:1, 0:1, type="n", xlab="", ylab="", axes=FALSE)
par(usr=c(0,1,0,1),xpd=NA)
legend(x="bottom",legend=c("PT","Nonpartisan","Anti-PT"),bty="n"
       ,pch=c(24,21,25),pt.bg=c(1,"white",1)
       ,col=1, xpd=NA,cex=1,ncol=3)

#############################################################
# Figure 1  
#############################################################

#Get the "correct" baselines for petista stereotypes
pt.NE <- 100*prop.table(svytable(~ptID+strNE, design=SurveyDesign),1)[2,2]
pt.CO <- 100*prop.table(svytable(~ptID+strCO, design=SurveyDesign),1)[2,2]
pt.02SM <- 100*prop.table(svytable(~ptID+str02SM, design=SurveyDesign),1)[2,2]
pt.0510SM <- 100*prop.table(svytable(~ptID+str0510SM, design=SurveyDesign),1)[2,2]
pt.atheist <- 100*prop.table(svytable(~ptID+stratheist, design=SurveyDesign),1)[2,2]
pt.catholic <- 100*prop.table(svytable(~ptID+strcatholic, design=SurveyDesign),1)[2,2]
pt.women <- 100*prop.table(svytable(~ptID+strwomen, design=SurveyDesign),1)[2,2]
pt.black <- 100*prop.table(svytable(~ptID+strblack, design=SurveyDesign),1)[2,2]

# Get the "correct" baselines for antipetista stereotypes
antipt.SSE <- 100*prop.table(svytable(~antiptID+strSSE, design=SurveyDesign),1)[2,2]
antipt.20SM <- 100*prop.table(svytable(~antiptID+str20SM, design=SurveyDesign),1)[2,2]
antipt.protestant <- 100*prop.table(svytable(~antiptID+strprotestant
                                             , design=SurveyDesign),1)[2,2]
antipt.catholic <- 100*prop.table(svytable(~antiptID+strcatholic, design=SurveyDesign),1)[2,2]
antipt.men <- 100*prop.table(svytable(~antiptID+strmen, design=SurveyDesign),1)[2,2]
antipt.white <- 100*prop.table(svytable(~antiptID+strwhite, design=SurveyDesign),1)[2,2]
antipt.over60 <- 100*prop.table(svytable(~antiptID+strover60 , design=SurveyDesign),1)[2,2]
antipt.3035 <- 100*prop.table(svytable(~antiptID+str3035 , design=SurveyDesign),1)[2,2]

truecategories <- ls()[grep("pt\\.",ls())]
truevalues <- sapply(truecategories,get)
print(as.matrix(truevalues))


## so that we can compute *actual* differences
antipt.NE <- 100*prop.table(svytable(~antiptID+strNE, design=SurveyDesign),1)[2,2]
antipt.CO <- 100*prop.table(svytable(~antiptID+strCO, design=SurveyDesign),1)[2,2]
antipt.02SM <- 100*prop.table(svytable(~antiptID+str02SM, design=SurveyDesign),1)[2,2]
antipt.0510SM <- 100*prop.table(svytable(~antiptID+str0510SM, design=SurveyDesign),1)[2,2]
antipt.atheist <- 100*prop.table(svytable(~antiptID+stratheist, design=SurveyDesign),1)[2,2]
antipt.women <- 100*prop.table(svytable(~antiptID+strwomen, design=SurveyDesign),1)[2,2]
antipt.black <- 100*prop.table(svytable(~antiptID+strblack, design=SurveyDesign),1)[2,2]

# Get the "correct" baselines for petista stereotypes
pt.SSE <- 100*prop.table(svytable(~ptID+strSSE, design=SurveyDesign),1)[2,2]
pt.20SM <- 100*prop.table(svytable(~ptID+str20SM, design=SurveyDesign),1)[2,2]
pt.protestant <- 100*prop.table(svytable(~ptID+strprotestant
                                         , design=SurveyDesign),1)[2,2]
pt.men <- 100*prop.table(svytable(~ptID+strmen, design=SurveyDesign),1)[2,2]
pt.white <- 100*prop.table(svytable(~ptID+strwhite, design=SurveyDesign),1)[2,2]
pt.over60 <- 100*prop.table(svytable(~ptID+strover60 , design=SurveyDesign),1)[2,2]
pt.3035 <- 100*prop.table(svytable(~ptID+str3035 , design=SurveyDesign),1)[2,2]

true.ptcategories <- ls()[grep("^pt\\.",ls())]
true.antiptcategories <- ls()[grep("^antipt\\.",ls())]
true.differences <- sapply(true.ptcategories,get)-sapply(true.antiptcategories,get)
names(true.differences) <- gsub("^pt\\.","",names(true.differences))
print(as.matrix(true.differences ))


# Just stereotypical categories
true.ptcategories <- ls()[grep("^pt\\.",ls())][c(1,6,10,15,5,11,3,12,14,9)]
true.antiptcategories <- ls()[grep("^antipt\\.",ls())][c(1,6,10,15,5,11,3,12,14,9)]
true.differences <- sapply(true.ptcategories,get)-sapply(true.antiptcategories,get)
names(true.differences) <- gsub("^pt\\.","",names(true.differences))
lab.true.differences <- c("Low Inc.","Non-white","NE","Women","Irreligious",
                          "Older","High Inc.","Protestant","White","Men")
pt.ordering <- order(true.differences,decreasing=F)
par(mar=c(5,6,1,1),mfrow=c(1,1))
plot(c(-25,25),c(1,length(true.differences))
     ,type="n",yaxt="n",bty="n"
     ,ylab="",xlab="",ylim=c(0.5,length(true.differences)+.5)
     ,main="True differences between petistas and antipetistas"
)
axis(side=2,at=1:length(true.differences)
     ,labels=lab.true.differences[pt.ordering]
     ,las=2,tick=F)
axis(side=1,at=c(-12.5,12.5)
     ,line=2,tick=F
     ,labels=c("Overrepresented in antipetistas","Overrepresented in petistas")
     ,cex=0.8)
abline(h=1:length(true.differences),col=gray(0.7),lty=3)
points(true.differences[pt.ordering],1:length(true.differences))
abline(v=0,col=gray(0.5))



#############################################################
# Figure 2 
#############################################################

# Using only stereotypical categories (dropping filler categories)
pt.truevalues <- truevalues[grep("^pt\\.",truecategories)][c(1,4,7,8,3)]
lab.pt.truevalues <- c("Low Inc.","Non-white","NE","Women","Irreligious")
antipt.truevalues <- truevalues[grep("^antipt\\.",truecategories)][c(5,1,6,8,4)]
lab.antipt.truevalues <- c("Older","High Inc.","Protestant","White","Men")

perceived.ptvals <- data.frame(rbind(
  extract.estse("profilespt_twosal"),
  extract.estse("profilespt_black"),
  extract.estse("profilespt_northeast"),
  extract.estse("profilespt_women"),
  extract.estse("profilespt_atheist")))

perceived.antiptvals <- data.frame(rbind(
  extract.estse("profilesantipt_older65"),
  extract.estse("profilesantipt_twentysal"),
  extract.estse("profilesantipt_protestant"),
  extract.estse("profilesantipt_white"),
  extract.estse("profilesantipt_men")))

par(mfrow = c(1,2),
    oma = c(2,0,0,0) + 0.1,
    mar = c(3,6,1,1) + 0.1)
pt.ordering <- order(perceived.ptvals$est-pt.truevalues)
plot(c(0,100),c(1,length(pt.truevalues))
     ,type="n",yaxt="n",bty="n",ylab="",xlab="",ylim=c(0.5,length(pt.truevalues)+.5))
axis(side=2,at=1:length(pt.truevalues)
     ,labels=lab.pt.truevalues[pt.ordering]
     ,las=2,tick=F)
abline(h=1:length(pt.truevalues),col=gray(0.7),lty=3)
points(pt.truevalues[pt.ordering],1:length(pt.truevalues))
segments(x0=perceived.ptvals$lb[pt.ordering],
         x1=perceived.ptvals$ub[pt.ordering],
         y0=1:length(pt.truevalues),
         y1=1:length(pt.truevalues))
points(perceived.ptvals$est[pt.ordering],1:length(pt.truevalues)
       ,pch=19)
mtext("Petista Stereotypes",line=2.2,side=1)

antipt.ordering <- order(perceived.antiptvals$est-antipt.truevalues)
par(mar=c(3,6,1,1))
plot(c(0,100),c(1,length(antipt.truevalues))
     ,type="n",yaxt="n",bty="n",ylab="",xlab="",ylim=c(0.5,length(antipt.truevalues)+.5))
axis(side=2,at=1:length(antipt.truevalues)
     ,labels=lab.antipt.truevalues[antipt.ordering]
     ,las=2,tick=F)
abline(h=1:length(antipt.truevalues),col=gray(0.7),lty=3)
points(antipt.truevalues[antipt.ordering],1:length(antipt.truevalues))
segments(x0=perceived.antiptvals$lb[antipt.ordering],
         x1=perceived.antiptvals$ub[antipt.ordering],
         y0=1:length(antipt.truevalues),
         y1=1:length(antipt.truevalues))
points(perceived.antiptvals$est[antipt.ordering],1:length(antipt.truevalues)
       ,pch=19)
mtext("Antipetista Stereotypes",line=2.2,side=1)
par(mfrow=c(1, 1), oma=rep(0, 4), mar=rep(0, 4), new=TRUE)
plot(0:1, 0:1, type="n", xlab="", ylab="", axes=FALSE)
legend(x="bottom",xpd=NA,legend=c("True value","Perceived value"),
       pch=c(1,19),lty=c(NA,1),bty="n",horiz=T)
rm(pt.ordering,antipt.ordering)

#################################################
# Figure 3 - True and perceived policy gaps
#################################################

issues <- c("envir","racial","militar",
            "condena","guns","inequality",
            "wedding",
            "gas","church","abortion")
issue.labels <- c("Environment","Racial justice",
                  "Limiting military appoint.",
                  "Disallow candidacies",
                  "Facilitate guns",
                  "Reduce inequality",
                  "Same sex marriage","Gas prices",
                  "Church tax exemptions","Abortion")
issue.invert <- c(T,T,T,F,F,T,T,T,F,T)

# Apply function that analyzes all issues
all.issues <- mapply(issue.sum,issues,issue.invert,issue.labels,SIMPLIFY=F)

real.differences <- t(sapply(all.issues,function(x){x$real.gap*100}))
print(real.differences)

perceived.differences <- t(sapply(all.issues,function(x){x$perceived.gap}))

ordering <- order(real.differences[,"est"])

par(mar=c(6.5,12,1,3),mfrow=c(1,1))
plot(c(-60,20),c(1,nrow(real.differences))
     ,type="n",yaxt="n",bty="n",ylab="",
     xlab="",ylim=c(0.5,nrow(real.differences)+.5),
     #main="Difference in support between petistas and antipetistas",
     xpd=NA)
axis(side=2,at=1:nrow(real.differences),
     labels=issue.labels[ordering],
     las=2,tick=F)
abline(h=1:nrow(real.differences),col=gray(0.7),lty=3)
abline(v=0,col=gray(0.7))
segments(x0=real.differences[ordering,"cilow"],
         x1=real.differences[ordering,"cihigh"],
         y0=1:nrow(real.differences),
         y1=1:nrow(real.differences))
points(real.differences[ordering,"est"],1:nrow(real.differences),bg="white",pch=21)
points(perceived.differences[ordering,"est"],1:nrow(perceived.differences)
       ,pch=19)
legend(x="bottom",xpd=NA,
       legend=c("True differences","Perceived differences"),
       pch=c(21,19),lty=c(1,NA),pt.bg=c("white",1),
       bty="n",horiz=T,inset=-0.25)


#################################################
# Figure D1 - Appendix 
#################################################

# Apply function that analyzes all issues
all.issues.pol <- mapply(issue.pol,issues,issue.invert,issue.labels,SIMPLIFY=F)

# Produce table
all.issues.pol <- do.call("rbind",all.issues.pol)

knitr::kable(all.issues.pol[order(all.issues.pol[,"pol"],decreasing=T),]
             ,caption = "Issue Polarization")

tmp <- merge(all.issues.pol,real.differences, by=0)
par(mar=c(4,4,1,1))
plot(tmp$pol,-tmp$est,bty="n",
     xlab=c("Issue Polarization Index"),
     ylab=c("Diferences in support across groups"),type="n")
tmp.reg <- lm(-est~pol,data=tmp)
abline(reg=tmp.reg,col=gray(0.5))
text(tmp$pol,-tmp$est,labels=tmp$Row.names,cex=0.8)


##########################################################
# Create APB measure
########################################################

df_labels$apb_pt <-  1/5 * (
  df_labels$profilespt_northeast - pt.NE  +
    df_labels$profilespt_twosal - pt.02SM  +
    df_labels$profilespt_atheist - pt.atheist  +
    df_labels$profilespt_women - pt.women  +
    df_labels$profilespt_black -  pt.black )

df_labels$apb_antipt <-  1/5 * (
  df_labels$profilesantipt_twentysal -  antipt.20SM +
    df_labels$profilesantipt_white -  antipt.white  + 
    df_labels$profilesantipt_protestant -  antipt.protestant  +
    df_labels$profilesantipt_men - antipt.men  +
    df_labels$profilesantipt_older65  -  antipt.over60 )


#####################################################################
# Figure 4 (4a and 4b) - Average Perception Bias by Group Membership
#####################################################################

typeapb_pt <- lm_robust(apb_pt~type5-1,data=df_labels,weights=df_labels$`weights#all`)
typeapb_antipt <- lm_robust(apb_antipt~type5-1,data=df_labels,weights=df_labels$`weights#all`)

pred_pt <- predict(typeapb_pt,interval="confidence"
                   ,newdata=data.frame(type5=unique(df_labels$type5)[c(5,3,1,4,2)]))

pred_antipt <- predict(typeapb_antipt,interval="confidence"
                       ,newdata=data.frame(type5=unique(df_labels$type5)[c(5,3,1,4,2)]))

xlabels <- c("Strong\nPetista","Petista","Nonpartisan",
             "Antipetista","Strong\nAntipetista")
nn <- length(xlabels)
par(mfrow=c(1,1),mar=c(6,4,1,3))
plot(c(1,nn),c(0,40),xaxt="n",bty="n",type="n",xlab="",ylab=""
     ,xpd=NA)
abline(h=pred_pt$fit[3,1],lty=3)
plotCI(1:nn,pred_pt$fit[,"fit"],
       ui=pred_pt$fit[,"upr"],li=pred_pt$fit[,"lwr"]
       ,col=1,lwd=1.5
       ,pch=19,ylab="",xlab="",yaxt="n",xaxt="n",bty="n",add=T)
axis(side=1,at=1:nn,labels=xlabels ,las=2)


par(mfrow=c(1,1),mar=c(6,4,1,3))
plot(c(1,nn),c(0,40),xaxt="n",bty="n",type="n",xlab="",ylab=""
     #,main="Antipetista Stereotypes"
     ,xpd=NA)
abline(h=pred_antipt$fit[3,1],lty=3)
plotCI(1:nn,pred_antipt$fit[,"fit"],
       ui=pred_antipt$fit[,"upr"],li=pred_antipt$fit[,"lwr"]
       ,col=1,lwd=1.5
       ,pch=19,ylab="",xlab="",yaxt="n",xaxt="n",bty="n",add=T)
axis(side=1,at=1:nn,labels=xlabels ,las=2)	


##################################################################################
# Figure 5 (5a and 5b) - Average Perception Bias by Level of Interest in Politics
##################################################################################

df_labels$interest <- df_labels$intpoliticalpond=="Muito interessado"|
  df_labels$intpoliticalpond=="Interessado"
df_labels$interest.low <- df_labels$intpoliticalpond=="Nada interessado"

## Estimate
reg.interest <- lm_robust(interest  ~type7-1,data=df_labels,weights=df_labels$`weights#all`)
reg.interest.low <- lm_robust(interest.low  ~type7-1,data=df_labels,weights=df_labels$`weights#all`)

# By interest 
intapb_pt <- lm_robust(apb_pt~intpoliticalpond-1,data=df_labels,weights=df_labels$`weights#all`)
intapb_antipt <- lm_robust(apb_antipt~intpoliticalpond-1,data=df_labels,weights=df_labels$`weights#all`)

pred_pt <- predict(intapb_pt,interval="confidence"
                   ,newdata=data.frame(intpoliticalpond=levels(df_labels$intpoliticalpond)[c(4,3,2,1)]))
pred_antipt <- predict(intapb_antipt,interval="confidence"
                       ,newdata=data.frame(intpoliticalpond=levels(df_labels$intpoliticalpond)[c(4,3,2,1)]))

xlabels <- gsub("\\s","\n",levels(df_labels$intpoliticalpond)[c(4,3,2,1)])

xlabels_english<-c("No\nInterest", "Low\nInterest", "Some\nInterest", "High\nInterest")

par(mfrow=c(1,1),mar=c(6,4,1,3))
plot(c(1,4),c(0,40),xaxt="n",bty="n",type="n",xlab="",ylab=""
     # ,main="Petista Stereotypes"
     ,xpd=NA)
abline(h=pred_pt$fit[1,"fit"],lty=2)
plotCI(1:4,pred_pt$fit[,"fit"],
       ui=pred_pt$fit[,"upr"],
       li=pred_pt$fit[,"lwr"]
       ,col=1,lwd=1.5
       ,pch=19,ylab="",xlab="",yaxt="n",xaxt="n",bty="n",add=T)
axis(side=1,at=1:4,labels=xlabels_english ,las=2)	

par(mfrow=c(1,1),mar=c(6,4,1,3))
plot(c(1,4),c(0,40),xaxt="n",bty="n",type="n",xlab="",ylab=""
     # ,main="Antipetista Stereotypes"
     ,xpd=NA)
abline(h=pred_antipt$fit[1,"fit"],lty=2)
plotCI(1:4,pred_antipt$fit[,"fit"],
       ui=pred_antipt$fit[,"upr"],li=pred_antipt$fit[,"lwr"]
       ,col=1,lwd=1.5
       ,pch=19,ylab="",xlab="",yaxt="n",xaxt="n",bty="n",add=T)
axis(side=1,at=1:4,labels=xlabels_english ,las=2)	



#############################################################
# Analyses: Regression tables
#############################################################

# Compute and table number of answer to each issue
knitr::kable(apply(is.na(sapply(all.issues,function(x){x$recoded.pt}))==F,2,sum)
             , caption = "Number of answers to each question")

# Assemble a data.frame with the recoded answers to each policy
# Recoded is "estimated share of each group taking the conservative position on the issue"
# These are extracted from "all.issues", the output of the estimating function 
recoded.answers.pt <- sapply(all.issues,function(x){x$recoded.pt})
recoded.answers.antipt <- sapply(all.issues,function(x){x$recoded.antipt})
colnames(recoded.answers.pt)<-paste(colnames(recoded.answers.pt),"pt",sep=".")
colnames(recoded.answers.antipt)<-paste(colnames(recoded.answers.antipt),"antipt",sep=".")

# Create a dataset with all answers in long format with individual indicator
issues.set.pt <- data.frame(id=df_labels$idEntrevista,
                            type=df_labels$type7,
                            interest=df_labels$interest,
                            w=df_labels$`weights#all`,
                            apb_pt=df_labels$apb_pt,
                            p_conservative=recoded.answers.pt)

issues.set.antipt <- data.frame(id=df_labels$idEntrevista,
                                type=df_labels$type7,
                                interest=df_labels$interest,
                                w=df_labels$`weights#all`,
                                apb_antipt=df_labels$apb_antipt,
                                p_conservative=recoded.answers.antipt)

# Reshape to create dataset for stacked analysis 
library(reshape2)

issues.set.pt <- na.omit(
  melt(issues.set.pt,id.vars=c("id","type","interest","apb_pt","w")))

issues.set.pt$extreme.partisan <- issues.set.pt$type=="antipetista.extreme"|
  issues.set.pt$type=="petista.extreme"
issues.set.pt$antipetista <- issues.set.pt$type=="antipetista"|
  issues.set.pt$type=="antipetista.extreme"|
  issues.set.pt$type=="antipetista.strong"

issues.set.antipt <- na.omit(
  melt(issues.set.antipt,id.vars=c("id","type","interest","apb_antipt","w")))
issues.set.antipt$extreme.partisan <- issues.set.antipt$type=="antipetista.extreme"|
  issues.set.antipt$type=="petista.extreme"
issues.set.antipt$antipetista <- issues.set.antipt$type=="antipetista"|
  issues.set.antipt$type=="antipetista.extreme"|
  issues.set.antipt$type=="antipetista.strong"

# Create the basic type in the stacked sets
issues.set.antipt$type3 <- gsub("\\..*","",issues.set.antipt$type)
issues.set.pt$type3 <- gsub("\\..*","",issues.set.pt$type)

recoded.gap <- abs(recoded.answers.pt-recoded.answers.antipt)
colnames(recoded.gap)<-paste(colnames(recoded.answers.pt),"antip.v.pt",sep=".")

# Create a dataset with all answers in long format with individual indicator
# Differently than preceding analysis, here we create only one dataset as gaps included both pt and antipt
issues.set.gap <- data.frame(id=df_labels$idEntrevista,
                             type=df_labels$type7,
                             w=df_labels$`weights#all`,
                             apb_pt=df_labels$apb_pt,
                             apb_antipt=df_labels$apb_antipt,
                             interest=df_labels$interest,
                             p_conservative_gap=recoded.gap)

# Reshape to create dataset for stacked analysis 
library(reshape2)
issues.set.gap <- na.omit(
  reshape2::melt(issues.set.gap,id.vars=c("id","type","apb_pt","apb_antipt","interest","w")))

issues.set.gap$extreme.partisan <- issues.set.gap$type=="antipetista.extreme"|
  issues.set.gap$type=="petista.extreme"
issues.set.gap$antipetista <- issues.set.gap$type=="antipetista"|
  issues.set.gap$type=="antipetista.extreme"|
  issues.set.gap$type=="antipetista.strong"
issues.set.gap$petista <- issues.set.gap$type=="petista"|
  issues.set.gap$type=="petista.extreme"|
  issues.set.gap$type=="petista.strong"

# Create the basic type in the stacked sets (we didn't need this before)
# We can delete the extreme and antipetista variables,above
issues.set.gap$type3 <- gsub("\\..*","",issues.set.gap$type)
issues.set.gap$type3 <- gsub("\\..*","",issues.set.gap$type)

true.shares.pt <- t(sapply(all.issues,function(x){summary(x$reg.agree)$coef["typepetista",]}))
true.shares.antipt <- t(sapply(all.issues,function(x){summary(x$reg.agree)$coef["typeantipetista",]}))

#############################################################
# Table 3 - Partisans Perceive Greater Policy Polarization
#############################################################

fe.gap.nofe <- lm_robust(value~type3,
                         clusters=id,
                         data=issues.set.gap,
                         se_type = "stata",
                         weights=w)

fe.gap <- lm_robust(value~type3,
                    fixed_effects=~variable,
                    clusters=id,
                    data=issues.set.gap,
                    se_type = "stata",
                    weights=w)

fe.gap.c <- lm_robust(value~type3+interest,
                      fixed_effects=~variable,
                      clusters=id,
                      data=issues.set.gap)


print(texreg(list(fe.gap.nofe,fe.gap,fe.gap.c)
             , include.ci = FALSE
             , custom.model.names = c("(1)","(2)","(3)") 
             , custom.coef.map = list("(Intercept)"="(Intercept)",
                                      "type3petista"="Petista",
                                      "type3antipetista"="Antipetista",
                                      "interestTRUE"="High interest")
             , include.pvalues = T
             , digits = 2
             , stars = c(0.001, 0.01, 0.05, 0.1)
             , label = "tab:preceivedpolarization"
             , caption = "Partisans Perceive Greater Policy Polarization"
             , caption.above =T
             , bold = 0.1
             , ci.force = FALSE))


#############################################################
# Table 4 - People Who Stereotype More Perceive Greater...
#############################################################

fe.pt <- lm_robust(value~apb_pt,
                   fixed_effects=~variable,
                   clusters=id,
                   data=issues.set.gap,
                   se_type = "stata",
                   weights=w)

fe.pt.t <- lm_robust(value~apb_pt+type3,
                     fixed_effects=~variable,
                     clusters=id,
                     data=issues.set.gap)

fe.pt.c <- lm_robust(value~apb_pt+type3+interest,
                     fixed_effects=~variable,
                     clusters=id,
                     data=issues.set.gap)

fe.antipt <- lm_robust(value~apb_antipt,
                       fixed_effects=~variable,
                       clusters=id,
                       data=issues.set.gap,
                       se_type = "stata",
                       weights=w)

fe.antipt.t <- lm_robust(value~apb_antipt+type3,
                         fixed_effects=~variable,
                         clusters=id,
                         data=issues.set.gap,
                         se_type = "stata",
                         weights=w)

fe.antipt.c <- lm_robust(value~apb_antipt+type3+interest,
                         fixed_effects=~variable,
                         clusters=id,
                         data=issues.set.gap,
                         se_type = "stata",
                         weights=w)

fe.both <- lm_robust(value~apb_pt+apb_antipt,
                     fixed_effects=~variable,
                     clusters=id,
                     data=issues.set.gap,
                     se_type = "stata",
                     weights=w)

fe.both.t <- lm_robust(value~apb_pt+apb_antipt+type3,
                       fixed_effects=~variable,
                       clusters=id,
                       data=issues.set.gap,
                       se_type = "stata",
                       weights=w)

fe.both.c <- lm_robust(value~apb_pt+apb_antipt+type3+interest,
                       fixed_effects=~variable,
                       clusters=id,
                       data=issues.set.gap,
                       se_type = "stata",
                       weights=w)

fe.0 <- lm_robust(value~type3+interest,
                  fixed_effects=~variable,
                  clusters=id,
                  data=issues.set.gap,
                  se_type = "stata",
                  weights=w)

fe.gap.nofe <- lm_robust(value~type3,
                         #fixed_effects=~variable,
                         clusters=id,
                         data=issues.set.gap,
                         se_type = "stata",
                         weights=w)

fe.gap <- lm_robust(value~type3,
                    fixed_effects=~variable,
                    clusters=id,
                    data=issues.set.gap,
                    se_type = "stata",
                    weights=w)

fe.gap.c <- lm_robust(value~type3+interest,
                      fixed_effects=~variable,
                      clusters=id,
                      data=issues.set.gap)

print(texreg(list(fe.pt,fe.pt.t,fe.pt.c,
                  fe.antipt,fe.antipt.t,fe.antipt.c
                  #fe.both,fe.both.t,fe.both.c
)
, include.ci = FALSE
, custom.model.names = c("(1)","(2)","(3)","(4)","(5)","(6)") 
, custom.coef.map = list("apb_pt"="APB^{PT}",
                         "apb_antipt"="APB^{antiPT}",
                         "type3petista"="Petista",
                         "type3antipetista"="Antipetista",
                         "interestTRUE"="High interest")
, include.pvalues = T
, digits = 2
, stars = c(0.001, 0.01, 0.05, 0.1)
, label = "tab:abpdifferences"
, caption.above =T
, bold = 0.1
, ci.force = FALSE))


#############################################################
# Table 5 - People with Strongest Partisan Group Stereotypes
#############################################################

fe.pt <- lm_robust(value~apb_pt,
                   fixed_effects=~variable,
                   clusters=id,
                   data=issues.set.pt,
                   se_type = "stata",
                   weights=w)

fe.pt.t <- lm_robust(value~apb_pt+type3,
                     fixed_effects=~variable,
                     clusters=id,
                     data=issues.set.pt,
                     se_type = "stata",
                     weights=w)

fe.pt.c <- lm_robust(value~apb_pt+type3+interest,
                     fixed_effects=~variable,
                     clusters=id,
                     data=issues.set.pt)

fe.pt.0 <- lm_robust(value~type3,
                     fixed_effects=~variable,
                     clusters=id,
                     data=issues.set.pt)

fe.antipt <- lm_robust(value~apb_antipt,
                       fixed_effects=~variable,
                       clusters=id,
                       data=issues.set.antipt,
                       se_type = "stata",
                       weights=w)

fe.antipt.t <- lm_robust(value~apb_antipt+type3,
                         fixed_effects=~variable,
                         clusters=id,
                         data=issues.set.antipt,
                         se_type = "stata",
                         weights=w)

fe.antipt.c <- lm_robust(value~apb_antipt+type3+interest,
                         fixed_effects=~variable,
                         clusters=id,
                         data=issues.set.antipt,
                         se_type = "stata",
                         weights=w)

fe.antipt.0 <- lm_robust(value~type3,
                         fixed_effects=~variable,
                         clusters=id,
                         data=issues.set.antipt,
                         se_type = "stata",
                         weights=w)

print(texreg(list(fe.pt,fe.pt.t,fe.pt.c,fe.antipt,fe.antipt.t,fe.antipt.c)
             , include.ci = FALSE
             , custom.header = list("Proportion PT Conservative"=1:3,"Proportion AntiPT Conservative"=4:6)
             , custom.model.names = c("(1)","(2)","(3)","(4)","(5)","(6)") 
             , custom.coef.map = list("apb_pt"="APBPT composition",
                                      "apb_antipt"="APBantiPT composition",
                                      "type3petista"="Petista",
                                      "type3antipetista"="Antipetista",
                                      "interestTRUE"="High interest")
             , include.pvalues = T
             , digits = 2
             , stars = c(0.001, 0.01, 0.05, 0.1)
             , label = "tab:replictab5"
             , caption = "People with the Most Prototype-Biased Party Perceptions See Partisans as More Likely to Take Party-Consistent Policy Positions"
             , caption.above =T
             , bold = 0.1
             , ci.force = FALSE))


#####################################################################
# Table 6 - Stereotyping is associated with desire for greater... 
#####################################################################

# Recode a couple of variables 
df_labels$extreme.partisan <- df_labels$type7=="antipetista.extreme"|
  df_labels$type7=="petista.extreme"
df_labels$antipetista <-  df_labels$type=="antipetista"|
  df_labels$type=="antipetista.extreme"|
  df_labels$type=="antipetista.strong"

# Recode social distance questions to numeric  
df_labels <- df_labels %>% mutate(
  #distance to antipt items
  dist_to_antipt_weddingpo = dplyr::recode(distanceoption1_weddingpo, `Muito feliz`=-2,`Feliz`=-1,`Indiferente`=0,`Infeliz`=1,`Muito infeliz`=2,`Não sei`=as.numeric(NA), .default = as.numeric(NA_character_)), 
  dist_to_antipt_workpol = dplyr::recode(distanceoption1_workpol, `Muito feliz`=-2,`Feliz`=-1,`Indiferente`=0,`Infeliz`=1,`Muito infeliz`=2,`Não sei`=as.numeric(NA), .default = as.numeric(NA_character_)), 
  dist_to_antipt_neightborpol = dplyr::recode(distanceoption1_neightborpolSQ002, `Muito feliz`=-2,`Feliz`=-1,`Indiferente`=0,`Infeliz`=1,`Muito infeliz`=2,`Não sei`=as.numeric(NA), .default = as.numeric(NA_character_)), 
  dist_to_antipt_distinctionpol = dplyr::recode(distanceoption1_distinctionpol, `Muito feliz`=-2,`Feliz`=-1,`Indiferente`=0,`Infeliz`=1,`Muito infeliz`=2,`Não sei`=as.numeric(NA), .default = as.numeric(NA_character_)), 
  #Distance to PT items
  dist_to_pt_weddingpo = dplyr::recode(distanceoption2_weddingpo, `Muito feliz`=-2,`Feliz`=-1,`Indiferente`=0,`Infeliz`=1,`Muito infeliz`=2,`Não sei`=as.numeric(NA), .default = as.numeric(NA_character_)), 
  dist_to_pt_workpol = dplyr::recode(distanceoption2_workpol, `Muito feliz`=-2,`Feliz`=-1,`Indiferente`=0,`Infeliz`=1,`Muito infeliz`=2,`Não sei`=as.numeric(NA), .default = as.numeric(NA_character_)), 
  dist_to_pt_neightborpol = dplyr::recode(distanceoption2_neightborpolSQ002, `Muito feliz`=-2,`Feliz`=-1,`Indiferente`=0,`Infeliz`=1,`Muito infeliz`=2,`Não sei`=as.numeric(NA), .default = as.numeric(NA_character_)), 
  dist_to_pt_distinctionpol = dplyr::recode(distanceoption2_distinctionpol, `Muito feliz`=-2,`Feliz`=-1,`Indiferente`=0,`Infeliz`=1,`Muito infeliz`=2,`Não sei`=as.numeric(NA), .default = as.numeric(NA_character_)))#end rename

antipt.socdist.vars <- c("dist_to_antipt_weddingpo",#casamento com antipetista
                         "dist_to_antipt_workpol",#trabalhar com antipetista
                         "dist_to_antipt_neightborpol",# placa bolsonaro presidente
                         "dist_to_antipt_distinctionpol")# condecoracao

pt.socdist.vars <- c("dist_to_pt_weddingpo",#casamento com  petista
                     "dist_to_pt_workpol",#trabalhar com  petista
                     "dist_to_pt_neightborpol",# placa lula presidente
                     "dist_to_pt_distinctionpol")# condecoracao pra lula

# Compute average distance to antipetistas/bolsonaro 
tmp <- subset(df_labels,select=c(antipt.socdist.vars))
alpha.antipt <- cronbach.alpha(na.omit(tmp))
df_labels$socdist.to.antipt <- apply(tmp,1,mean)#original index without imputation

##Number of missing items (216 cases missing 1 item)
missing.items.antipt <- apply(is.na(tmp),1,sum,na.rm=T)
print(table(missing.items.antipt))

# Compute average distance to petistas/lula  
tmp <- subset(df_labels,select=c(pt.socdist.vars))
alpha.pt <- cronbach.alpha(na.omit(tmp))
df_labels$socdist.to.pt <- apply(tmp,1,mean)


## Looking only at petistas relative to antipetistas
df_petistas <- subset(df_labels,type=="petista")
reg.towards.antipetistas <- lm_robust(socdist.to.antipt~apb_antipt,
                                      data=df_petistas,
                                      weights=df_petistas$`weights#all`)
reg.towards.antipetistas.c <- lm_robust(socdist.to.antipt~apb_antipt+extreme.partisan+interest,
                                        data=df_petistas,
                                        weights=df_petistas$`weights#all`)

## Looking only at antipetstas relative to petistas
df_antipetistas <- subset(df_labels,type=="antipetista")
reg.towards.petistas <- lm_robust(socdist.to.pt~apb_pt,
                                  data=df_antipetistas,
                                  weights=df_antipetistas$`weights#all`)
reg.towards.petistas.c <- lm_robust(socdist.to.pt~apb_pt+extreme.partisan+interest,
                                    data=df_antipetistas,
                                    weights=df_antipetistas$`weights#all`)

# Pool the data and focus on "outgroup"
df_labels$socdist_outgroup <- ifelse(df_labels$type=="petista",df_labels$socdist.to.antipt,
                                     ifelse(df_labels$type=="antipetista",df_labels$socdist.to.pt,
                                            NA))
df_labels$abp_outgroup <- ifelse(df_labels$type=="petista",df_labels$apb_antipt,
                                 ifelse(df_labels$type=="antipetista",df_labels$apb_pt,
                                        NA))

regp.01 <- lm_robust(socdist_outgroup~abp_outgroup,
                     data=df_labels,
                     weights=df_labels$`weights#all`)
regp.01c <- lm_robust(socdist_outgroup~abp_outgroup+extreme.partisan+antipetista+interest,
                      data=df_labels,
                      weights=df_labels$`weights#all`)
regp.00 <- lm_robust(socdist_outgroup~extreme.partisan+antipetista+interest,
                     data=df_labels,
                     weights=df_labels$`weights#all`)

print(texreg(list(reg.towards.petistas,reg.towards.petistas.c,
                  reg.towards.antipetistas,reg.towards.antipetistas.c,regp.01,regp.01c)
             , include.ci = FALSE
             , custom.header = list("Rel to PT"=1:2,"Rel to AntiPT"=3:4,"Rel to outgroup"=5:6)
             , custom.model.names = c("(1)","(2)","(3)","(4)","(5)","(6)") 
             , custom.coef.map = list("apb_pt"="APB PT composition",
                                      "apb_antipt"="APB antiPT composition",
                                      "abp_outgroup"="APB outgroup composition",
                                      "extreme.partisanTRUE"="Extreme partisan",
                                      "antipetistaTRUE"="Antipetista",
                                      "interestTRUE"="High interest")
             , include.pvalues = T
             , digits = 3
             , stars = c(0.001, 0.01, 0.05, 0.1)
             , label = "tab:replictab6"
             , caption = "Individuals with greater APB are more socially distant from outgroup (positive values mean more distance)"
             , caption.above =T
             , bold = 0.1
             , ci.force = FALSE))

#############################################################
# Table D1 - Appendix
#############################################################

# Apply function that analyzes all issues
all.issues.pol <- mapply(issue.pol,issues,issue.invert,issue.labels,SIMPLIFY=F)

# Produce table
all.issues.pol <- do.call("rbind",all.issues.pol)

knitr::kable(all.issues.pol[order(all.issues.pol[,"pol"],decreasing=T),]
             ,caption = "Issue Polarization")
